#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Wed Jun 12 05:18:31 2019

@author: hyungmokson
"""

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Mon May 27 17:15:28 2019

@author: hyungmokson
"""
import matplotlib as mpl
from matplotlib import rcParams
from matplotlib import rc
import matplotlib.pyplot as plt

mpl.use("pgf")

#from matplotlib.backends.backend_pgf import FigureCanvasPgf
#mpl.backend_bases.register_backend('pdf', FigureCanvasPgf)
#
print("MATPLOTLIB STUFF:")
print("matplotlib version: " + mpl.__version__)
print("matplotlib directory: " + mpl.__file__)
print(" ")
def figsize(scale):
    fig_width_pt = 350.                          # Get this from LaTeX using \the\textwidth
    inches_per_pt = 1.0/72.27                       # Convert pt to inch
    golden_mean = 0.8#(np.sqrt(5.0)-1.0)/2.0            # Aesthetic ratio (you could change this)
    fig_width = fig_width_pt*inches_per_pt*scale    # width in inches
    fig_height = fig_width*golden_mean              # height in inches
    fig_size = [fig_width,fig_height]
    return fig_size

#
pgf_with_latex = {                      # setup matplotlib to use latex for output
    "pgf.texsystem": "pdflatex",        # change this if using xetex or lautex
    "text.usetex": True,                # use LaTeX to write all text
#    "font.family": "serif",
#    "font.serif": [],                   # blank entries should cause plots to inherit fonts from the document
#    "font.sans-serif": [],
#    "font.monospace": [],
    "axes.linewidth": 1,
    "axes.labelsize": 7,               # LaTeX default is 10pt font.
#    "text.fontsize": 10,
    "legend.fontsize": 7,               # Make the legend/label fonts a little smaller
    "xtick.labelsize": 7,
    "ytick.labelsize": 7,
    "figure.figsize": figsize(0.9),     # default fig size of 0.9 textwidth
    "pgf.preamble": [
#        r"\usepackage[utf8x]{inputenc}",    # use utf8 fonts becasue your computer can handle it :)
#        r"\usepackage[T1]{fontenc}",
        r"\usepackage{amsmath}",
        r"\usepackage{amstext}", 
        r"\usepackage[utf8x]{inputenc}",
        r"\usepackage[T1]{fontenc}",
        r"\usepackage{cmbright}",# plots will be generated using this preamble
        ]
    }

pgf_with_latex_simple = {                      # setup matplotlib to use latex for output
    "pgf.texsystem": "pdflatex",        # change this if using xetex or lautex
    "text.usetex": True,
#    "font.family": "serif",
#    "font.serif": [],                   # blank entries should cause plots to inherit fonts from the document
#    "font.sans-serif": [],
#    "font.monospace": [],
    "axes.labelsize": 7,               # LaTeX default is 10pt font.
#    "text.fontsize": 10,
    "legend.fontsize": 7,               # Make the legend/label fonts a little smaller
    "xtick.labelsize": 7,
    "ytick.labelsize": 7              # use LaTeX to write all text]
    }

rcParams.update(pgf_with_latex)

rcParams['font.family'] = 'sans-serif'

###########################################################################################################
###########################################################################################################
from matplotlib.ticker import FuncFormatter

import numpy as np
from numpy import linalg
import matplotlib.pyplot as plt
import os
from os import listdir
from os.path import isfile, join
from scipy.optimize import curve_fit
from scipy import interpolate
from scipy.interpolate import interp1d
from scipy.interpolate import CubicSpline

if __name__ == "__main__":

    working_dir = os.getcwd()
    sub = os.listdir(working_dir)
    print(working_dir)
    print(sub)
    
    main= join(working_dir, 'single evap')
    main2= join(working_dir, 'double evap')
      
    def outputData(main, fileName, flg):
        if flg == 0:
            path = join(main, 'without Na')
        elif flg == 1:
            path = join(main, 'with Na')
        else:
            path = join(main, 'base cond')
#            path = join(main, 'base cond_new')
        
        path = join(path, 'csv')
        name = str(fileName) + '.csv'
        temp = np.genfromtxt(join(path, name), delimiter=',', dtype=float, comments = "#", names=('var', 'data', 'err'))
        x = np.array(temp['var'])
        y = np.array(temp['data'])
        e = np.array(temp['err'])
                    
        output = []
        output.append(x)
        output.append(y)
        output.append(e)
        return output

    plotpath = join(working_dir, 'plot save') 
             
    def trapdepthToPower(trapdepth):
        return 0.3063 * trapdepth - 0.0172
    
    print("")
    print("trap power for 0.25V = " + str(trapdepthToPower(0.25)) + " Watt")
    
    ## axial trap freq for Na is 33.8 KHz for 4 Volt
    ## The trap freq ratio between molecule/Na = 1.4387, given the polarizability ratio of 2.6
    def trapfreq_mol_ax(trapdepth):
        f_ax_4volt = 33800 * 1.4387
        return f_ax_4volt * np.sqrt(trapdepthToPower(trapdepth)/trapdepthToPower(4.))
    
    def trapfreq_na_ax(trapdepth):
        f_ax_4volt = 33800
        return f_ax_4volt * np.sqrt(trapdepthToPower(trapdepth)/trapdepthToPower(4.))
    
    def trapfreq_mol(trapdepth):
        return np.sqrt((800. *650.)/trapdepthToPower(4.5) * trapdepthToPower(trapdepth))

    def trapdepth_uk_mol(trapdepth):
        return 0.25 * (29 *10**-3/(6.02*1.38)) * (50*10**-6)**2 * trapfreq_mol(trapdepth)**2 *1e6 * (2*3.141592)**2 

    def trapfreq_na(trapdepth):
        return np.sqrt((600. * 460.)/trapdepthToPower(5.) * trapdepthToPower(trapdepth))

    def trapdepth_uk_na(trapdepth):
        return 0.25 * (23 *10**-3/(6.02*1.38)) * (50*10**-6)**2 *trapfreq_na(trapdepth)**2 *1e6 * (2*3.141592)**2

#    print(trapfreq_mol(0.5))
    setpt_list = np.array([0.25, 0.4, 0.5, 0.62, 0.75, 1, 1.5, 2, 4, 5])
    print("trap power = " + str(trapdepthToPower(setpt_list)) + " Watt")
    print("trap depth for MOL = " + str(trapdepth_uk_mol(setpt_list)) + " uK")
    print("for set point of " + str(setpt_list) + " Volt")
    
    print("")
    print(" ---- FOR MOLECULE ---- ")
    print("5 V trapdepth molecule = " + str(trapdepth_uk_mol(5)) + " uk")
    print("4 V trapdepth molecule = " + str(trapdepth_uk_mol(4)) + " uk")
    print("0.5 V trapdepth molecule = " + str(trapdepth_uk_mol(0.5)) + " uk")
    print("0.4 V trapdepth molecule = " + str(trapdepth_uk_mol(0.4)) + " uk")
    print("0.25 V trapdepth molecule = " + str(trapdepth_uk_mol(0.25)) + " uk")

    print("")
    print(" ---- FOR NA ---- ")
    print("5 V trapdepth Na = " + str(trapdepth_uk_na(5)) + " uk")
    print("0.5 V trapdepth na = " + str(trapdepth_uk_na(0.5)) + " uk")
    print("0.4 V trapdepth na = " + str(trapdepth_uk_na(0.4)) + " uk")
    print("0.25 V trapdepth na = " + str(trapdepth_uk_na(0.25)) + " uk")
    print("")
    
    a = (1596*1e-3)/2.
    sigma = 720. ## um
    l_eff = np.sqrt(8./np.pi) * sigma
    particle_per_site = l_eff/a
    
    plank_const = 6.62607*10**(-34)
    boltzman_const = 1.38065 * 10**(-23)
    num_factor = (plank_const/(np.sqrt(2)*boltzman_const))**3
    factor = trapfreq_mol(5.)**2 * trapfreq_mol_ax(5.) * num_factor
    factor = factor/particle_per_site
    new_factor = factor * (10**-6)**(-3) ## this is to take the "micro-Kelvin" unit
    print("")
    print("Calculated factor = " + str(new_factor))
    print("")
        
    factor = (1.23888e-3)/particle_per_site
    print("OLD factor = " + str(factor))
    print("")
    factor = new_factor
    print(factor)

    def calc_psd(num, temp, num_err, temp_err):
        psd = factor * num * temp**-3
        psd_err = factor * np.sqrt((num_err * temp**-3)**2 + (3 * num * temp_err * temp**-4)**2)        
        return psd, psd_err
    
    ###########################################################################
    ########## base cond ###############
    ###########################################################################    
        
    molnum_base = outputData(main, 'mol num', -1)    
    t_molnum_base = trapdepthToPower(molnum_base[0])
    n_molnum_base = molnum_base[1]
    err_molnum_base = molnum_base[2]
    
    moltemp_base = outputData(main, 'mol temp', -1)    
    t_moltemp_base = trapdepthToPower(moltemp_base[0])
    T_moltemp_base = moltemp_base[1]
    err_moltemp_base = moltemp_base[2]

    t_molpsd_base = t_moltemp_base
    psd_base, err_psd_base = calc_psd(n_molnum_base, T_moltemp_base, err_molnum_base, err_moltemp_base)
    molpsd_base = [t_molpsd_base, psd_base, err_psd_base]
    
    print("")
    print("Base mol num = " + str(n_molnum_base) + " +/- " + str(err_molnum_base))
    print("Base mol temp = " + str(T_moltemp_base) + " +/- " + str(err_moltemp_base))
    print("")
    
    psd_true_base = psd_base
    err_psd_true_base = err_psd_base
    
    baseflg = False
    
    ###########################################################################
    ########## For single evap ###############
    ###########################################################################
    molnum = outputData(main, 'mol num', 1)
    t_molnum = trapdepthToPower(molnum[0])
    n_molnum = molnum[1]
    err_molnum = molnum[2]
    
    molnum[0] = t_molnum
    
    moltemp = outputData(main, 'mol temp', 1)
    trapdepth = trapdepthToPower(moltemp[0])
    t_moltemp_avg = trapdepthToPower(moltemp[0])
    T_moltemp_avg = moltemp[1]
    err_moltemp_avg = moltemp[2] 

    moltemp[0] = t_moltemp_avg

    psd, err = calc_psd(n_molnum, T_moltemp_avg, err_molnum, err_moltemp_avg)
    molpsd = [trapdepth, psd, err]    
    
    molnum_no = outputData(main, 'mol num', 0)
    t_molnum_no = trapdepthToPower(molnum_no[0])
    n_molnum_no = molnum_no[1]
    err_molnum_no = molnum_no[2]
    
    molnum_no[0] = t_molnum_no
    
    moltemp_no = outputData(main, 'mol temp', 0)
    trapdepth_no = trapdepthToPower(moltemp_no[0])
    t_moltemp_avg_no = trapdepthToPower(moltemp_no[0])
    T_moltemp_avg_no = moltemp_no[1]
    err_moltemp_avg_no = moltemp_no[2] 

    moltemp_no[0] = t_moltemp_avg_no

    psd_no, err_no = calc_psd(n_molnum_no, T_moltemp_avg_no, err_molnum_no, err_moltemp_avg_no)
    molpsd_no = [trapdepth_no, psd_no, err_no]
        
    psd_base = []
    psd_base.append(trapdepthToPower(np.array(5)))
    psd_base.append(np.array(1127))
    psd_base.append(np.array(107))
    
    psd_true = psd
    err_psd_true = err
    
    psd_true_no = psd_no
    err_psd_true_no = err_no

    #######################################
    ##### Na num only for single evap #####
    #######################################
    nanum = outputData(main, 'na num', 1)
    t_nanum = nanum[0]
    n_nanum = nanum[1]
    err_nanum = nanum[2]    
    
    ###########################################################################
    ########## For double evap ###############
    ###########################################################################
    main = main2
    molnum2 = outputData(main, 'mol num', 1)
    t_molnum2 = trapdepthToPower(molnum2[0])
    n_molnum2 = molnum2[1]
    err_molnum2 = molnum2[2]
    
    molnum2[0] = t_molnum2
    
    moltemp2 = outputData(main, 'mol temp', 1)
    trapdepth2 = trapdepthToPower(moltemp2[0])
    t_moltemp_avg2 = trapdepthToPower(moltemp2[0])
    T_moltemp_avg2 = moltemp2[1]    
    err_moltemp_avg2 = moltemp2[2] 
    ## update for plotting
    moltemp2[0] = t_moltemp_avg2

    psd2, err2 = calc_psd(n_molnum2,T_moltemp_avg2, err_molnum2, err_moltemp_avg2)
    molpsd2 = [trapdepth2, psd2, err2]
    
    print("")
    print("trap depth in Volt = " + str(molpsd2[0]))
    trapdepth_mol_temporary = trapdepth_uk_mol(molpsd2[0])
    print("trap depth in uK for mol = " + str(trapdepth_mol_temporary) + " uk")
    print("trap depth in Watt = " + str(trapdepth2))
    print("")
    
    ######### True PSD of molecules for 5V lattice ##############
#    trapdepth_y_true, psd_y_true, err_y_true = outputData(main, 'psd_y_true')
      
    fig, ax = plt.subplots(3, 1, sharex=True, facecolor="white")
    ax[0].tick_params(axis='both',  which='both')
    ax[1].tick_params(axis='both',  which='major')
    ax[2].tick_params(axis='both',  which='major')
    
    hor = 8.9/2.54
    fig.set_size_inches(hor, 9.5/6. * hor)
    fig.subplots_adjust(hspace=0.05)
    
    lw = 1.
    ms = 3.8
    lw_list = [1., 1., 1.]
    fmt_list = ['o', 's', 'x']
    facecolor_list = ['white','#ff0000', '#000000']
    color_list=['#0000ff','#ff0000', '#000000']

    toPlot = []
    toPlot.append(molpsd2)
    toPlot.append(molpsd)
    toPlot.append(molpsd_no)
    if baseflg == True:
        toPlot.append(molpsd_base)
        
    annotate_pos = (-60, 170)
#    toPlot.append(psd_base)

    ax_index_num = 0
    ax_index_temp = 1        
    ax_index_psd = 2
    i = 0
    for data in toPlot:
        x = data[0]
        y = data[1]
        err = data[2]
        
        (_, caps, _) = ax[ax_index_psd].errorbar(x,y, yerr = err, fmt = fmt_list[i], mew=lw_list[i], capsize=0, markersize=ms, c = color_list[i], markerfacecolor = facecolor_list[i]) ##
        i += 1
    for cap in caps:
        cap.set_markeredgewidth(0.)

    ax[ax_index_psd].set_ylabel("Phase space density")
    ax[ax_index_psd].set_ylim([0.8e-4, 3e-2])
    ax[ax_index_psd].set_yscale('log')
    
    def funcForEyes(amp_factor, x, y, p):
        x_new = x
        y_new = y * amp_factor
        cs = CubicSpline(x_new, y_new, bc_type=['not-a-knot', 'natural'], extrapolate = True)
        interp_value = cs(p)
        return interp_value

    guide_a = 0.15
#    guide_lw = 1.65
    guide_lw = 1.55

    i = 0
    interp_original = False
    indexes = []
    ## for double evap
    index1 = np.ones(len(toPlot[0][0]))
    if interp_original is False:
        index1[0] = 0.95
#        index1[1] = 1.1
        index1[2] = 1.1
        index1[3] = 1.
        index1[-2] = 1.02
    indexes.append(index1)
    
    ## for single evap
    index2 = np.ones(len(toPlot[1][0]))
    if interp_original is False:
        index2[1] = 1.1
        index2[5] = 1.05
        index2[6] = 1.05
        index2[7] = 1.
        index2[8] = 1.5
    indexes.append(index2)

    ## for NO evap
    index3 = np.ones(len(toPlot[2][0]))
    if interp_original is False:
        index3[0] = 0.9
        index3[1] = 1.05
        index3[3] = 0.9
        index3[4] = 1.1
        index3[5] = 0.85
        index3[6] = 1.05
        index3[7] = 0.9
        index3[-1] = 1.1
    indexes.append(index3)
    for data in toPlot:
        x = data[0]
        y = data[1]
        err = data[2]
        
        p = np.arange(np.min(toPlot[i][0]), np.max(toPlot[i][0]), 0.001)
        guide_line = funcForEyes(indexes[i], x, y, p)
    
        ax[ax_index_psd].plot(p, guide_line, color=color_list[i], alpha=guide_a, linewidth=guide_lw)
        i += 1
 
    ## draw dotted-line for the initial psd
    x = np.linspace(0, 10., 10)
    y = np.ones(10) * psd_true_base 
    lower_lim = psd_true_base - molpsd_base[2]
    upper_lim = psd_true_base + molpsd_base[2]
    y1 = np.ones(10) * lower_lim
    y2 = np.ones(10) * upper_lim
#    ax[0].plot(x, y1, 'k--')
#    ax[0].plot(x, y2, 'k--')
    ax[ax_index_psd].plot(x, y, 'k--', linewidth = 0.65)
    ax[ax_index_psd].axhspan(lower_lim, upper_lim, xmin = 0, xmax = 10, alpha = 0.15, facecolor = 'k')

    ####################
    ##### mol num ######
    ####################
    toPlot = []
    toPlot.append(molnum2)
    toPlot.append(molnum)
    toPlot.append(molnum_no)
    if baseflg == True:
        toPlot.append(molnum_base)
    
    num_factor = 1e4
    i = 0
    interp_original = False
    indexes = []
    ## for double evap
    index1 = np.ones(len(toPlot[0][0]))
    if interp_original is False:
        index1[0] = 1.05
        index1[1] = 1.01
        index1[2] = 1.05
        index1[3] = 1.02
    indexes.append(index1)
    
    ## for single evap
    index2 = np.ones(len(toPlot[1][0]))
    if interp_original is False:
        index2[2] = 1.015
        index2[4] = 1.01
        index2[5] = 0.855
        index2[6] = 1.005
        index2[7] = 1.005
    indexes.append(index2)

    ## for NO evap
    index3 = np.ones(len(toPlot[2][0]))
    if interp_original is False:
        index3[0] = 0.98
        index3[1] = 0.95
        index3[2] = 1.02
        index3[5] = 1.05
        index3[6] = 1.05
        index3[7] = 0.85
        index3[8] = 1.2
    indexes.append(index3)

    i = 0
    for data in toPlot:
        x = data[0]
        y = data[1]/num_factor
        err = data[2]/num_factor
        
        (_, caps, _) = ax[ax_index_num].errorbar(x,y, yerr = err, fmt = fmt_list[i], mew=lw_list[i], capsize=0, markersize=ms, c = color_list[i], markerfacecolor = facecolor_list[i]) ##
        i += 1

    for cap in caps:
        cap.set_markeredgewidth(0.)
    
    legend_box = ['Double evap.', 'Single evap.', 'Single evap.\n (without Na)']    
    if baseflg == True:
        legend_box.append('No evap.')
        
    ax[0].legend(legend_box, frameon = True, ncol = 3, facecolor = 'w', handlelength = 0.0, columnspacing = 1., borderpad = 0.1, edgecolor ='w',  fontsize = 7, framealpha = 0., markerscale = 1)
#loc='upper left',
#labelspacing = 0.55
    ax[ax_index_num].set_ylabel("Molecule num. (" + str(r"$\times 10^{4}$") + ")")

    i = 0
    for data in toPlot:
        x = data[0]
        y = data[1]/num_factor
        err = data[2]/num_factor

        p = np.arange(np.min(toPlot[i][0]), np.max(toPlot[i][0]), 0.001)
        guide_line = funcForEyes(indexes[i], x, y, p)
    
        ax[ax_index_num].plot(p, guide_line, color=color_list[i], alpha=guide_a, linewidth=guide_lw)
#        ax[1].annotate(r"$ $", xy=(340, 95), xycoords='axes points', size=11, ha='right', va='top')

        i += 1
    
    ####################
    ##### mol temp ######
    ####################
    toPlot = []
    toPlot.append(moltemp2)
    toPlot.append(moltemp)
    toPlot.append(moltemp_no)
    if baseflg == True:
        toPlot.append(moltemp_base)
    i = 0
    interp_original = False
    indexes = []
    ## for double evap
    index1 = np.ones(len(toPlot[0][0]))
    if interp_original is False:
        index1[0] = 0.95

    indexes.append(index1)
    
    ## for single evap
    index2 = np.ones(len(toPlot[1][0]))
    if interp_original is False:
        index2[1] = 0.97
        index2[5] = 0.9
        index2[6] = 1.0
        index2[7] = 1.05
        index2[8] = 0.98
    indexes.append(index2)

    ## for NO evap
    index3 = np.ones(len(toPlot[2][0]))
    if interp_original is False:
        index3[0] = 1.02
        index3[1] = 0.98
        index3[4] = 0.935
        index3[5] = 1.07

    indexes.append(index3)

    for data in toPlot:
        x = data[0]
        y = data[1]
        err = data[2]
        
        (_, caps, _) = ax[ax_index_temp].errorbar(x,y, yerr = err, fmt = fmt_list[i], mew=lw_list[i], capsize=0, markersize=ms, c = color_list[i], markerfacecolor = facecolor_list[i]) ##        
#        ax[ax_index_temp].annotate(r'$\textbf{b.}', xy= (-60, 167), xycoords='axes points', size=14, ha='left', va='top')
#        ax[ax_index_temp].annotate(r"$ $", xy=(365, 95), xycoords='axes points', size=14, ha='right', va='top')

        p = np.arange(np.min(toPlot[i][0]), np.max(toPlot[i][0]), 0.001)
        guide_line = funcForEyes(indexes[i], x, y, p)
    
        ax[ax_index_temp].plot(p, guide_line, color=color_list[i], alpha=guide_a, linewidth=guide_lw)

        i += 1
    for cap in caps:
    #cap.set_color('red')
        cap.set_markeredgewidth(0.)
    
    ylabel = "Molecule temp. " + str(r"$(\mu\mathrm{K})$")
    ax[ax_index_temp].set_ylabel(ylabel)
    
    ## draw dotted-line for the initial psd
    
    x = np.linspace(0, 10., 10)
    y = np.ones(10) * T_moltemp_base
    lower_lim = T_moltemp_base - err_moltemp_base
    upper_lim = T_moltemp_base + err_moltemp_base
    y1 = np.ones(10) * lower_lim
    y2 = np.ones(10) * upper_lim
#    ax[ax_index_temp].plot(x, y1, 'k--')
#    ax[ax_index_temp].plot(x, y2, 'k--')
    ax[ax_index_temp].plot(x, y, 'k--', linewidth = 0.65)
    ax[ax_index_temp].axhspan(lower_lim, upper_lim, xmin = 0, xmax = 10, alpha = 0.15, facecolor = 'k')
    
    for i in [0,1,2]:
        ax[i].yaxis.set_ticks_position('both')
        ax[i].xaxis.set_ticks_position('both')
        ax[i].tick_params(which='both', direction='in')  
#        ax[i].grid(which='both', linestyle='dotted')
        
        ax[i].set_xscale('log')    
#        ax[i].set_xlim([0.45, 5.5])
        ax[i].set_xlim([0.085, 2.3])
        if i == ax_index_temp:
#            ax[i].set_xticks([0.4, 0.5, 0.7, 1, 2, 3, 5])
            ax[i].set_ylim([0.35, 3.55])
            ax[i].set_yticks([0.5, 1.5, 2.5, 3.5])

        if i == ax_index_psd:
            ax[i].set_xticks([0.1, 0.2, 0.3, 0.4, 0.5, 1., 2.])
            ax[i].get_xaxis().set_major_formatter(mpl.ticker.ScalarFormatter())
            ax[i].set_xlabel("1596 nm laser power (Watt)")
           
        if i == ax_index_num:
            ax[i].set_ylim([0.15, 1.65])
            ax[i].set_yticks([0.25, 0.65, 1.05, 1.45])
                
    s = 8
    fig.canvas.draw()
    labels = [r'$\textbf{a.}', r'$\textbf{b.}', r'$\textbf{c.}']
    axlabels = [fig.text(0,0, label, fontsize=s, fontweight="bold", va="top", ha="left") for a, label in zip(ax, labels)]

    def update_labels(evt=None):
        trans = fig.transFigure.inverted()
        i = 0
        for a, label in zip(ax, axlabels):
            bbox = a.get_tightbbox(fig.canvas.get_renderer())
            if i == 0:
                x = bbox.x0
#            x = bbox.x0
            y = bbox.y1
            label.set_position(trans.transform_point([x, y]))
            i += 1
            
    update_labels()

    if baseflg == True:
        plt.savefig(join(plotpath, "psd_avg_single_double_with_base.pdf")) 
    else:
        plt.savefig(join(plotpath, "psd_avg_single_double.pdf"), bbox_inches = 'tight', pad_inches=0.01, dpi = 1000)
        
    fig.show()    
    i_max1 = np.where(molpsd[1]==np.max(molpsd[1]))[0][0]
    i_max2 = np.where(molpsd2[1]==np.max(molpsd2[1]))[0][0]    
    print("")
    print("the initial PSD = " + str(psd_true_base) + " / " + str(err_psd_true_base))
    print("the max PSD for SINGLE evap  = " + str(np.max(molpsd[1])) + " / " + str(molpsd[2][i_max1]))
    print("the max PSD for double evap = " + str(np.max(molpsd2[1])) + " / " + str(molpsd2[2][i_max2]))
    print("")   
    maxsingle_div_single = np.max(molpsd[1])/psd_true[-1]
    err = maxsingle_div_single * np.sqrt((molpsd[2][i_max1]/np.max(molpsd[1]))**2 + (err_psd_true[-1]/psd_true[-1])**2)
    print("max psd/5V-psd for single evap = " + str(maxsingle_div_single) + " / " + str(err))
    maxsingle_div_noevap =  np.max(molpsd[1])/psd_true_base
    err2 = maxsingle_div_noevap * np.sqrt((molpsd[2][i_max1]/np.max(molpsd[1]))**2 + (err_psd_true_base/psd_true_base)**2)
    print("max psd/no-evap-psd for single evap = " + str(maxsingle_div_noevap) + " / " + str(err2))
    print("")
    maxdouble_div_double = np.max(molpsd2[1])/psd_true[-1]
    err3 = maxdouble_div_double * np.sqrt((molpsd[2][i_max2]/np.max(molpsd2[1]))**2 + (err_psd_true[-1]/psd_true[-1])**2)
    print("max psd/5V-psd for double evap = " + str(maxdouble_div_double) + " / " + str(err3))
    maxdouble_div_noevap = np.max(molpsd2[1])/psd_true_base
    err4 = maxdouble_div_noevap * np.sqrt((molpsd2[2][i_max2]/np.max(molpsd2[1]))**2 + (err_psd_true_base/psd_true_base)**2)
    print("max psd/no-evap-psd for double evap = " + str(maxdouble_div_noevap) + " / " + str(err4))
#    plt.show()
